In [1]:
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPoint
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt

import folium
from folium.plugins import MarkerCluster
import geopandas as gpd

import numpy as np

import gurobipy as gp
from gurobipy import GRB

Lab 2: Fire Station Service Coverage Analysis¶

Overview¶

In this lab, we will explore the concept of spatial efficiency in the context of facility placement. Spatial efficiency refers to the optimal allocation of resources in geographical space to maximize coverage while minimizing redundancy. This concept is crucial in various fields such as urban planning, logistics, and environmental management.

The primary objective of this lab is to analyze the spatial efficiency of existing and optimized facility placement strategies. We will investigate how different factors, such as coverage distance and the number of stations, influence spatial efficiency metrics.

Objectives¶

  • Understand the concept of spatial efficiency in facility placement.
  • Analyze the spatial distribution of existing facilities and their coverage areas.
  • Optimize facility placement to improve spatial efficiency metrics.
  • Evaluate the impact of coverage distance and the number of stations on spatial efficiency.

Tools and Technologies¶

  • Python programming language
  • Geospatial libraries (e.g., GeoPandas, Folium)
  • Gurobipy for spatial optimization
  • Jupyter Lab environment

Dataset¶

  • Santa Barbara County Faces: polygon shapefile with Land / Water Flag
  • Santa Barbara County Fire Stations: point shapefile

Lab Structure¶

  1. Data Exploration and Preprocessing
  2. Analysis of Existing Facility Placement
  3. Optimization of Facility Placement
  4. Evaluation of Spatial Efficiency Metrics
  5. Final Report: Conclusion and Discussion

Step 1. Data Exploration and Preprocessing¶

In [2]:
# In this cell, let's read data

# Data source: granted by SB Fire District GIS Manager
stations = gpd.read_file("Data/Fire_Stations_SBC.shp")

# Data source: https://catalog.data.gov/dataset/tiger-line-shapefile-2019-county-santa-barbara-county-ca-topological-faces-polygons-with-all-ge
county = gpd.read_file("tl_2019_06111_faces/tl_2019_06111_faces.shp")
In [3]:
county.plot("LWFLAG")
Out[3]:
<Axes: >
No description has been provided for this image
In [4]:
# Unfortunately, county data contains the ocean part, which we don't want..
    # check it out in QGIS or ArcGIS if you want!

# Let's remove the water part by using the column "LWFLAG". 
    # When LWFLAG column value is L, it means the geometry is part of land,
    # While LWFLAG column value is W, it means the geometry is waterbody
county = county.loc[county['LWFLAG'] == 'L'].reset_index(drop=True)
county.plot("LWFLAG")
Out[4]:
<Axes: >
No description has been provided for this image
In [ ]:
 
In [5]:
# Let's unary_union every geometry in county GeoDataFrame
    # Because it contains all the small blocks, but we want the whole county boundary.

county_shape = county.unary_union
    # Note: unary_union returns a GEOMETRY only!
county_shape
Out[5]:
No description has been provided for this image
In [6]:
county_shape.geom_type
Out[6]:
'MultiPolygon'
In [7]:
# Plot the two datasets in one plot
# Does it look good..?

A = stations.plot(color="red")
# To plot the GEOMETRY, county_shape with another gpd.GeoDataFrame, stations, 
# let's convert it into gpd.GeoSeries
A = gpd.GeoSeries(county_shape).plot()
stations.plot(ax=A)

# When the plot doesn't look plausible, first thing to do is to look into the CRS.
Out[7]:
<Axes: >
No description has been provided for this image
No description has been provided for this image
In [8]:
# Let's check stations dataset's crs
    # It's in a Projected Coordiante system, where the unit is "foot"
stations.crs
Out[8]:
<Projected CRS: EPSG:2229>
Name: NAD83 / California zone 5 (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura.
- bounds: (-121.42, 32.76, -114.12, 35.81)
Coordinate Operation:
- name: SPCS83 California zone 5 (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
In [9]:
# Let's check county dataset's crs
    # It's in a Geographic Coordinate system, where the unit is "degree"
county.crs #This is taking the angular unit. Makes it complicated. So we should match the 
#county and stations crs.
Out[9]:
<Geographic 2D CRS: EPSG:4269>
Name: NAD83
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands.
- bounds: (167.65, 14.92, -40.73, 86.45)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
In [10]:
# Let's convert the county dataset to stations crs. 
# Because stations have a planar coordinate system (foot linear unit).
county_ = county.to_crs(stations.crs)
# Actually, what this line does is to create a new GeoDataFrame
    # with one geometry, unary_union of water county blocks, in stations.crs.
county_shape = gpd.GeoDataFrame(geometry=[county.to_crs(stations.crs).unary_union], crs=stations.crs)
county_shape
Out[10]:
geometry
0 MULTIPOLYGON (((6118330.643 1541434.591, 61182...
In [11]:
# Now, let's plot them together.

A = stations.plot(color="red", zorder=100, markersize=10)
county_shape.plot(ax=A, facecolor='None', edgecolor="purple")

# Review for SCALEBAR!
scalebar = ScaleBar(dx=1, scale_formatter=lambda data, unit: f'{data},000 ft', length_fraction=0.25)
A.add_artist(scalebar)
Out[11]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x178b66d90>
No description has been provided for this image
In [12]:
county_shape.set_crs(county_.crs)
Out[12]:
geometry
0 MULTIPOLYGON (((6118330.643 1541434.591, 61182...
In [22]:
# How about interactive mapping?

import folium
from folium.plugins import MousePosition
import geopandas as gpd


county_shape_4326 = county_shape.to_crs(4326)

# Calculate the centroid of the polygon
centroid = county_shape_4326.geometry.centroid.values[0].coords[0]

# Create a Folium map centered at the centroid of the polygon
m = folium.Map(location=[centroid[1], centroid[0]], zoom_start=9)

# Add the polygon to the map
folium.GeoJson(county_shape_4326, style_function=lambda x: { 'fillColor': 'transparent'}).add_to(m)

folium.GeoJson(stations.to_crs(4326)).add_to(m)

# Add MousePosition control to display cursor coordinates
#### Use the cursor coordinates to define your Area of Interest! ####
MousePosition().add_to(m)

# Display the map
m
/var/folders/0f/r37lnqs969n4hggs14l43w280000gn/T/ipykernel_42483/3721739741.py:11: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroid = county_shape_4326.geometry.centroid.values[0].coords[0]
Out[22]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [23]:
# Following 4 float values are given for lab question 1.
# For lab question 2, you need to change it to your own boundary coordinates
max_x, min_x = -118.933753, -118.669013
max_y, min_y = 34.316016, 34.169861

#34.206016, -118.923753
#34.135861, -118.769013

# Create the rectangle geometry, which forms are study region boundary.
boundary = Polygon([(min_x, min_y), (min_x, max_y), (max_x, max_y), (max_x, min_y), (min_x, min_y)])

# Add the boundary to the map
folium.GeoJson(boundary.__geo_interface__, style_function=lambda x: {'color': 'red', 'fillColor': 'transparent'}).add_to(m)

m
Out[23]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [ ]:
 
In [24]:
# Do we agree with the study region boundary?
# Let's check the area

# let's have it as a GeoSeries, with crs defined (4326)
# And convert it back to feet crs (stations.crs)
boundary = gpd.GeoSeries(boundary, crs=4326).to_crs(stations.crs)[0]

# And retrieve the area value
boundary.area / 2.788e+7
# According to google, let's divide the value with 2.788e+7 to convert square ft to square miles
# The boundary is about 97 square miles
Out[24]:
152.64019759003588
In [25]:
#boundary[0]

Step 2. Analysis of Existing Facility Placement¶

In [26]:
# To take a subset of fire stations within our study region, 
    # Let's get the boolean (True/False) array indicating whether each station is within our boundary.

stations.within(boundary).sum()
#.sum()
Out[26]:
15
In [27]:
# Remember? We can use .loc[] to take the subset!
    # Then we use reset_index(drop=True) to make the row numbers consecutive.
study_stations = stations.loc[stations.within(boundary)].reset_index(drop=True)

# Check out how many stations we have in our study region.
len(study_stations)
Out[27]:
15
In [28]:
county_shape
Out[28]:
geometry
0 MULTIPOLYGON (((6118330.643 1541434.591, 61182...
In [29]:
# Also, let's take an intersection of county polygon with our study region.

boundary_county = county_shape.intersection(boundary)
boundary_county
Out[29]:
0    MULTIPOLYGON (((6279557.580 1921454.930, 62797...
dtype: geometry
In [30]:
# Let's define how much fire station can cover!
# I would try 1 mile firstly
COVERAGE_DISTANCE = 1 * 5280 # because 1*5280 feet is 1 mile 

# Draw buffer from each fire station in study_stations GeoDataFrame 
    # Buffer distance is COVERAGE_DISTANCE we've defined
# stations_coverage = study_stations.buffer(COVERAGE_DISTANCE) 
stations_coverage = study_stations.buffer(COVERAGE_DISTANCE)

# Let's plot the county polgyon within our study region, and
    # stations_coverage buffers
A = stations_coverage.plot(facecolor="None")
boundary_county.plot(ax=A, zorder=-1)
plt.axis("off")

# How much area is covered?
    # To get to know, firstly take unary_union of the total station_coverage
total_stations_coverage = stations_coverage.unary_union
    # Then take the intersection with boundary_county, to exclude water area
total_stations_coverage = total_stations_coverage.intersection(boundary_county)
    # Then divide the area of the updated total_station_coverage with the boundary_county.area
coverage = total_stations_coverage.area / boundary_county.area
    # Multiply 100 to get the % value
coverage = coverage * 100

# Check out how much percentage of study region is covered by current fire stations.
    # Do they agree with each other? the coverage % and plot?
print(f'{round(coverage.values[0],2)} % covered')
26.37 % covered
No description has been provided for this image
In [31]:
total_stations_boundary(0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[31], line 1
----> 1 total_stations_boundary(0)

NameError: name 'total_stations_boundary' is not defined
In [ ]:
 
In [32]:
# This is a short example, how to create 2 pair combinations in python.

people = ["Jiwon", "Zhongqi", "Rosemary", "Rita"]

for i in range(len(people)-1):
    for j in range(i+1, len(people)):
        print((i, j), people[i],",", people[j])
        
# Just be familiarized with this syntax :)
    # The combination result makes sense, right?
(0, 1) Jiwon , Zhongqi
(0, 2) Jiwon , Rosemary
(0, 3) Jiwon , Rita
(1, 2) Zhongqi , Rosemary
(1, 3) Zhongqi , Rita
(2, 3) Rosemary , Rita
In [33]:
stations_coverage.reset_index(drop=True)
Out[33]:
0     POLYGON ((6299243.133 1890829.251, 6299217.708...
1     POLYGON ((6315022.851 1885664.677, 6314997.426...
2     POLYGON ((6304560.049 1902991.506, 6304534.624...
3     POLYGON ((6286542.009 1892333.313, 6286516.584...
4     POLYGON ((6336773.659 1889174.404, 6336748.234...
5     POLYGON ((6319012.314 1896007.759, 6318986.890...
6     POLYGON ((6295460.572 1920553.311, 6295435.147...
7     POLYGON ((6344539.690 1921791.230, 6344514.265...
8     POLYGON ((6301732.777 1927484.823, 6301707.352...
9     POLYGON ((6362490.000 1921820.000, 6362464.575...
10    POLYGON ((6318544.579 1914554.090, 6318519.154...
11    POLYGON ((6352715.852 1930322.205, 6352690.428...
12    POLYGON ((6336583.341 1928323.266, 6336557.916...
13    POLYGON ((6330451.435 1921355.349, 6330426.011...
14    POLYGON ((6285789.870 1891454.840, 6285764.445...
dtype: geometry
In [34]:
# How much redundant coverage?

redundant_coverage = 0
    # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
for i in range(len(stations_coverage)-1):
    # Take one station coverage buffer
    cover_i = stations_coverage[i]
    
    for j in range(i+1, len(stations_coverage)):
        # Take another station coverage buffer which hasn't paired with i
        cover_j = stations_coverage[j]
        print("\t combination: ", i, j)
        
        # With the two stations combination, calculate the intersection area.
            # and sum up!
        redundant_coverage += cover_i.intersection(cover_j).area

print() # To add blank line
print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
	 combination:  0 1
	 combination:  0 2
	 combination:  0 3
	 combination:  0 4
	 combination:  0 5
	 combination:  0 6
	 combination:  0 7
	 combination:  0 8
	 combination:  0 9
	 combination:  0 10
	 combination:  0 11
	 combination:  0 12
	 combination:  0 13
	 combination:  0 14
	 combination:  1 2
	 combination:  1 3
	 combination:  1 4
	 combination:  1 5
	 combination:  1 6
	 combination:  1 7
	 combination:  1 8
	 combination:  1 9
	 combination:  1 10
	 combination:  1 11
	 combination:  1 12
	 combination:  1 13
	 combination:  1 14
	 combination:  2 3
	 combination:  2 4
	 combination:  2 5
	 combination:  2 6
	 combination:  2 7
	 combination:  2 8
	 combination:  2 9
	 combination:  2 10
	 combination:  2 11
	 combination:  2 12
	 combination:  2 13
	 combination:  2 14
	 combination:  3 4
	 combination:  3 5
	 combination:  3 6
	 combination:  3 7
	 combination:  3 8
	 combination:  3 9
	 combination:  3 10
	 combination:  3 11
	 combination:  3 12
	 combination:  3 13
	 combination:  3 14
	 combination:  4 5
	 combination:  4 6
	 combination:  4 7
	 combination:  4 8
	 combination:  4 9
	 combination:  4 10
	 combination:  4 11
	 combination:  4 12
	 combination:  4 13
	 combination:  4 14
	 combination:  5 6
	 combination:  5 7
	 combination:  5 8
	 combination:  5 9
	 combination:  5 10
	 combination:  5 11
	 combination:  5 12
	 combination:  5 13
	 combination:  5 14
	 combination:  6 7
	 combination:  6 8
	 combination:  6 9
	 combination:  6 10
	 combination:  6 11
	 combination:  6 12
	 combination:  6 13
	 combination:  6 14
	 combination:  7 8
	 combination:  7 9
	 combination:  7 10
	 combination:  7 11
	 combination:  7 12
	 combination:  7 13
	 combination:  7 14
	 combination:  8 9
	 combination:  8 10
	 combination:  8 11
	 combination:  8 12
	 combination:  8 13
	 combination:  8 14
	 combination:  9 10
	 combination:  9 11
	 combination:  9 12
	 combination:  9 13
	 combination:  9 14
	 combination:  10 11
	 combination:  10 12
	 combination:  10 13
	 combination:  10 14
	 combination:  11 12
	 combination:  11 13
	 combination:  11 14
	 combination:  12 13
	 combination:  12 14
	 combination:  13 14

Redundant Coverage:  3.01068997223571  square miles
In [38]:
# Let's create a function to calculate a redundant coverage!

def calculate_redundant_coverage(coverage_buffers):
    redundant_coverage = 0
        # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
    for i in range(len(coverage_buffers)-1):
        # Take one station coverage buffer
        cover_i = coverage_buffers[i]

        for j in range(i+1, len(coverage_buffers)):
            # Take another station coverage buffer which hasn't paired with i
            cover_j = coverage_buffers[j]

            # With the two stations combination, calculate the intersection area.
                # and sum up!
            redundant_coverage += cover_i.intersection(cover_j).area

    print() # To add blank line
    print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
    return redundant_coverage / 2.788e+7
In [39]:
current_rc = calculate_redundant_coverage(stations_coverage)
Redundant Coverage:  3.01068997223571  square miles

Step 3. Optimization of Facility Placement¶

In [40]:
# Check this GeoDataFrame attribute, bounds!
    # It's convenient to find minx, miny, maxx, maxy
boundary.bounds
Out[40]:
(6279223.634048559, 1884871.2206143597, 6359655.31843054, 1938695.8015625787)
In [41]:
# Let's save them as variables.

minx, miny, maxx, maxy = boundary.bounds
print(minx, miny, maxx, maxy)
6279223.634048559 1884871.2206143597 6359655.31843054 1938695.8015625787
In [51]:
# Here, we are going to create fishnet points.
    # Fishnet points are points spaced with an equal interval.
    # You'll see when you see the output!

# Define the interval between points 
#### Note: Adjust interval value as needed for lab question 2
interval = 5280/4  # default: 5280/5 feet (1 mile / 5 = 0.2 mile) 


# Create arrays of x and y coordinates using np.arange
x_coords = np.arange(minx, maxx, interval)
y_coords = np.arange(miny, maxy, interval)

# Create a list to store the points
fishnet_points = []

# Generate points for the fishnet
for y in y_coords:
    for x in x_coords:
        fishnet_points.append((x, y))

# Print the number of points generated
print("Number of points in the fishnet:", len(fishnet_points))

fishnet_points = gpd.GeoSeries([Point(pt_cd) for pt_cd in fishnet_points])
fishnet_points.plot(figsize=(15,5))
Number of points in the fishnet: 2501
Out[51]:
<Axes: >
No description has been provided for this image
In [45]:
fishnet_points
Out[45]:
0       POINT (6279223.634 1884871.221)
1       POINT (6280983.634 1884871.221)
2       POINT (6282743.634 1884871.221)
3       POINT (6284503.634 1884871.221)
4       POINT (6286263.634 1884871.221)
                     ...               
1421    POINT (6351383.634 1937671.221)
1422    POINT (6353143.634 1937671.221)
1423    POINT (6354903.634 1937671.221)
1424    POINT (6356663.634 1937671.221)
1425    POINT (6358423.634 1937671.221)
Length: 1426, dtype: geometry
In [46]:
boundary
Out[46]:
No description has been provided for this image
In [52]:
### Note: Here, we want the fishnet points only if 
    # they are within our study region, and within the existing stations coverage
    
#fishnet_points = fishnet_points.loc[fishnet_points.within intersects(boundary_county[0])]
fishnet_points = fishnet_points.loc[fishnet_points.intersects(total_stations_coverage[0])].reset_index(drop=True)

fishnet_points.plot(figsize=(15,5))
Out[52]:
<Axes: >
No description has been provided for this image
In [53]:
# Overlay the fishnet_points with existing fire stations

A = fishnet_points.plot(figsize=(15,5), color="black", markersize=10)
study_stations.plot(ax=A, color="blue", marker="*", markersize=400)
Out[53]:
<Axes: >
No description has been provided for this image

Location Set Covering Problem¶

Given:

  • $ n $ = number of demand points ($ \text{num\_demands} $)
  • $ m $ = number of facilities ($ \text{num\_facilities} $)
  • $ D_{ij} $ = distance between demand point $ i $ and facility $ j $
  • $ N_i $ = set of facilities that can cover demand point $ i $
  • $ x_j $ = binary decision variable indicating whether facility $ j $ is selected ($ 1 $) or not ($ 0 $)

Objective: $ \text{Minimize} \quad \sum_{j=1}^{m} x_j $

Subject to: $ \sum_{j \in N_i} x_j \geq 1 $ for $ i = 1, 2, \ldots, n $

Where:

  • $ \sum_{j \in N_i} $ denotes the sum over facilities in the set $ N_i $, i.e., the facilities that can cover demand point $ i $.
In [129]:
fruit_pair = {
    "Jiwon":"Pineapple",
    "Rita":"Orange",
    "Julie":"Mango",
}
In [130]:
fruit_pair["Julie"]
Out[130]:
'Mango'
In [55]:
#dict(zip(range(num_demands), [[] for _ in range(num_demands)]))
In [133]:
range(len(fishnet_points))
Out[133]:
range(0, 0)
In [135]:
# study_stations
In [56]:
# Can we cover the same amount with Fewer stations?
# Let's solve this optimization problem, LSCP !!

num_demands = len(fishnet_points) 
num_facilities = len(study_stations)

D_ij = [[np.sqrt((fishnet_points.x[i] - study_stations.geometry.x[j])**2 + (fishnet_points.y[i] - study_stations.geometry.y[j])**2) 
         for j in range(num_facilities)] for i in range(num_demands)]

N_i = dict(zip(range(num_demands), [[] for _ in range(num_demands)]))

for i in range(num_demands):
    for j in range(num_facilities):
        if D_ij[i][j] <= COVERAGE_DISTANCE:
            N_i[i].append(j)

# Create a new model
model = gp.Model("LSCP")

# Decision variables
# x[i] = 1 if facility i is selected, 0 otherwise
x = model.addVars(num_facilities, vtype=GRB.BINARY, name="x")

# Objective function: minimize the number of selected facilities
model.setObjective(gp.quicksum(x[j] for j in range(num_facilities)), GRB.MINIMIZE)

# Constraints
# Each demand point must be covered by at least one selected facility
for i in range(num_demands):
    model.addConstr(gp.quicksum(x[j] for j in N_i[i]) >= 1, name=f"cover_demand_point_{j}")

# Optimize the model
model.optimize()
Restricted license - for non-production use only - expires 2025-11-24
Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (mac64[arm] - Darwin 22.1.0 22A400)

CPU model: Apple M1
Thread count: 8 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 624 rows, 15 columns and 655 nonzeros
Model fingerprint: 0xa6167d99
Variable types: 0 continuous, 15 integer (15 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 15.0000000
Presolve removed 624 rows and 15 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)

Solution count 1: 15 

Optimal solution found (tolerance 1.00e-04)
Best objective 1.500000000000e+01, best bound 1.500000000000e+01, gap 0.0000%
In [136]:
N_i
Out[136]:
{}
In [57]:
selected_bool = [x[i].x > 0.5 for i in range(num_facilities)]
selected_facilities = study_stations.loc[selected_bool].reset_index(drop=True)
unselected_facilities = study_stations.loc[~np.array(selected_bool)].reset_index(drop=True)

print("Objective value: ", model.objVal, sum(selected_bool))
print("Current running stations: ", len(study_stations))
Objective value:  15.0 15
Current running stations:  15
In [58]:
A = gpd.GeoSeries(MultiPoint(fishnet_points)).plot(figsize=(15,5), color="grey", markersize=5)
unselected_facilities.plot(ax=A, color="blue", legend=True, label="Unselected", marker="*", markersize=300, zorder=10)
selected_facilities.plot(ax=A, color="red", marker="*", markersize=300, legend=True, label="Selected")
selected_facilities.buffer(COVERAGE_DISTANCE).plot(ax=A, edgecolor="red", facecolor="None", linestyle="dashed")
plt.legend()
plt.axis("off")
/var/folders/0f/r37lnqs969n4hggs14l43w280000gn/T/ipykernel_42483/901057089.py:2: UserWarning: The GeoDataFrame you are attempting to plot is empty. Nothing has been displayed.
  unselected_facilities.plot(ax=A, color="blue", legend=True, label="Unselected", marker="*", markersize=300, zorder=10)
Out[58]:
(6270866.863114576, 6366853.006545841, 1877623.8010640834, 1938363.0812950088)
No description has been provided for this image
In [59]:
# Convert GeoPandas GeoSeries to GeoDataFrame
fishnet_gdf = gpd.GeoDataFrame(geometry=fishnet_points, crs=stations.crs).to_crs(4326)

# Create a base map
m = folium.Map(location=[fishnet_gdf.unary_union.centroid.y, fishnet_gdf.unary_union.centroid.x], zoom_start=11)



# Plot fishnet
for idx, row in fishnet_gdf.iterrows():
    folium.CircleMarker(location=[row.geometry.y, row.geometry.x],
                        radius=1,
                        color='grey', fillcolor="grey").add_to(m)

# Plot study stations
for idx, row in unselected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='blue'),
                   popup=row['Label']).add_to(m)

# Plot selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='red', icon='star'),
                   popup=row['Label']).add_to(m)

# Plot buffer around selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.vector_layers.Circle(location=[row.geometry.y, row.geometry.x],
                                 radius=COVERAGE_DISTANCE/3.281, # feet to meter conversion
                                 color='red',
                                 fill=False,
                                 dash_array='10').add_to(m)

# Display the map
m
Out[59]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [60]:
# Calculate Redundant Coverage of selected facilities
    # Hint: you can use the function we defined!
    
optimized_rc = calculate_redundant_coverage(selected_facilities.buffer(COVERAGE_DISTANCE))
    
# Is redundant coverage decreased with optimization?
Redundant Coverage:  3.01068997223571  square miles

Step 4. Evaluation of Spatial Efficiency Metrics¶

In [61]:
# Let's create the spatial efficienty

# Minimum (ideal) number of facilities needed to cover the demands covered by current facilities
# divided by the number of current facilities (fire stations).

spatial_efficiency = model.objVal / len(study_stations) * 100

print("Spatial Efficienty is ", round(spatial_efficiency, 2), " %")

# Higher efficienty, closer to the ideal number
# Lower efficienty, more redundancy, further from the ideal number
Spatial Efficienty is  100.0  %
In [62]:
print("Coverage Distance (mi): " , COVERAGE_DISTANCE) # Double Check this value before putting into table ** 
print("Current Stations #: " ,len(study_stations))
print("Optimized Stations #: ", len(selected_facilities))
print("Current Redundant Coverage (sqmi): ", current_rc)
print("Optimized Redundant Coverage (sqmi): ", optimized_rc)
print("Spatial Efficiency (%): ", spatial_efficiency)
Coverage Distance (mi):  5280
Current Stations #:  15
Optimized Stations #:  15
Current Redundant Coverage (sqmi):  3.01068997223571
Optimized Redundant Coverage (sqmi):  3.01068997223571
Spatial Efficiency (%):  100.0

Final Report¶

  1. Fill out the table below.
    You can create a table by coding (pd.DataFrame) or in any other software (Exel, Google Docs ... etc). ( 8 pts )
Coverage Distance (mi) Current Stations # Optimized Stations # Current Redundant Coverage (sqmi) Optimized Redundant Coverage (sqmi) Spatial Efficiency (%)
1
2
3
4
  1. Discuss the findings from the Table. How can this be reflected into Urban/Environmental Systems? ( 3 pts )

  2. Start a new notebook. This notebook will be saved as an HTML page.
    When this quarter concludes, upload the HTML to your personal portfolio website.
    Your want your portfolio different from others.
    Therefore, choose your own study area (boundary) using the folium map with cursor location.
    Hint: Ensure that only an urban region is selected.

    Experiment with different COVERAGE_DISTANCE values as demonstrated in Question 1.
    You're encouraged to customize the provided code skeleton by copying, pasting, editing, or functionizing it to create your unique Fire Service Coverage Analysis Report!
    Submit your exported HTML as file. ( 7 points )

Be sure to keep a copy of your Lab 2 HTML and all related materials handy, just in case you want to make edits for your personal portfolio later on!!

You can export your notebook in HTML by File > Download as > HTML (.html) in your jupyter lab notebook.
In case it doesn't work.., we will figure out week 2 for the lab 2!

In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: